Stressed backbone and elasticity of random central-force systems 
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We use a new algorithm to find the stress-carrying backbone of "generic" site-diluted triangular 
lattices of up to 10® sites. Generic lattices can be made by randomly displacing the sites of a regular 
lattice (see Fig. 1). The percolation threshold is Pc ~ 0.6975±0.0003, the correlation length exponent 
v = 1.16 ±0.03 and the fractal dimension of the backbone Db = 1.78 ±0.02. The number of "critical 
bonds" (if you remove them rigidity is lost) on the backbone scales as , with x = 0.85 ±0.05. The 
Young's modulus is also calculated. 
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The forces between atoms can often be divided into two 
classes "central forces" and "angular forces" (e.g. cova- 
lent bonds). In engineering, structures composed of bars 
connected at nodes (e.g. some bridges), get their rigid- 
ity primarily from the tensile and compressive stiffness 
of the bars (these are central- force terms). Structures 
of this sort arc called "trusses", while those in which 
the angle forces (or beam-bending) are important are 
called "frames". It is simple to see that systems which 
are dominated by angle forces support an applied stress 
as long as they are simply connected. In contrast, sys- 
tems with only central forces require higher order con- 
nectivity, the simplest rigid structure being a triangle. 
In many applications; for example in granular media^, 
glasses^, gels^'* and in engineering design, the disorder in 
a central-force structure is important and must be con- 
sidered. The stress-bearing paths of central-force systems 
have been primarily studied by brute-force solution of the 
force equations^"^. Although useful and important, this 
method is slow and subject to roundoff errors for large 
structures. An efficient method for relating the connec- 
tivity of a central-force structure to its ability to carry 
stress is an important and, in general unsolved, problem. 
One exception to this is two-dimensional random lattices, 
for which exact conditions^""^^ relating connectivity to 
"rigidity" have existed for over a decade. However, till 
recently^^ there has been no efficient implementation of 
these conditions, and their associated algorithms in ei- 
ther physics or engineering. This paper and the preced- 
ing paper by Jacobs and Thorpe (JT)^^ describe the first 
implementations of these ideas. We use our algorithm, 
to calculate the stressed backbone, and in combination 
with an iterative solver, to find the elastic properties of 
these backbones. Wc also identify the critical (red) bonds 
as those whose removal would lead to loss of rigidity, and 
study their scaling properties. For reasons outlined below 
these methods apply to randomly displaced (or "generic" 
- see Fig. la) central- force lattices. 




FIG. 1. A configuration that is unstable to shear on a reg- 
ular lattice, but is stable on a displaced lattice (dotted lines in- 
dicate absent bonds), a) The configuration in the "bar-joint" 
representation (28 joints and 53 bars), b) The configuration 
in the "body-bar" representation (2 bodies and 3 bars). 

Our ability to determine, from connectivity informa- 
tion alone, whether a central-force structure contains a 
stress carrying path is based on Laman's theorem^. 
Laman's theorem A random lattice (see below for a 
precise definition) consisting of N nodes and B bonds so 



that 2N — B = 3 is rigid if and only if there is no subset 
of the lattice, consisting of n nodes connected by b bonds, 
for which 2n — 6 < 3 is violated. 

This is the "bar-joint" statement of Laman's theorem. 
The origin of the expression 2n — 6 = 3 is easy to under- 
stand. Each node (joint) in two dimensions has two de- 
grees of freedom (two translations) , and each bond (bar) 
is a constraint (for example in Fig. la, n = 28, 6 = 53). 
In the expression 2n — 6 = 3, the 3 is there because in two 
dimensions a rigid body (in this case the whole lattice or 
cluster) has 3 degrees of freedom (two translations and 
a rotation). 2n — 6 = 3 is the two dimensional version 
of a general constraint counting argument introduced by 



Maxwell. 
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FIG. 2. The rigidity threshold as a function of sample 
size. AS with periodic boundary conditions (*), AS with open 
boundary conditions(-l-), IS with open boundary conditions 
(x) IS with periodic boundary conditions(o). The lattice sizes 

(L) (number of configurations) used are as follows; 16(2 x 10^), 
32(10^), 64(8 X 10"*), 128(2 x 10"*), 256(1.2 x 10''), 512(2 x 10^), 
1024(2 X 10^). 

However the new feature here is that constraint count- 
ing is exact in two dimensions provided it is implemented 
at all length scales (Unfortunately this result does not 
extend to three dimsensions, where counterexamples^^ 
to the three-dimensional extension of this argument, 
3n — 6 = 6, are known to exist). However, even in 2- 
d a naive algorithm must check all subclusters of a set 
of N nodes and so is not polynomial complete. How- 
ever Laman's theorem may be implemented by using 
the "bipartite matching" algorithm from graph theory 
which, when refined as described below, scales as N^'^^ 
for finding the stressed backbone at the rigidity percola- 
tion point. 

Our implementation of Laman's theorem is a cluster la- 
belling algorithm^ ^. Although we do site percolation, 
where p is the probability that a site is occupied, the 
algorithm works by testing a newly added bond against 
the configuration of rigid clusters already on the lattice. 
For a given p, we find the site configuration, and from it 
the configuration of present bonds. Then wc start with 
an empty lattice and add the present bonds one at a time. 
Each rigid cluster is a "body" with 3 degrees of freedom, 
so we must generalise the statement 2n — 6 = 3 of the 
original lattice to 3n(,od — 6 = 3, where njyod is the num- 



ber of bodies (or rigid clusters) in a configuration. For 
example the configuration of Fig. la, has two bodies and 
3 bars (see Fig. lb). A key component of the algorithm 
is the realisation by Hendrickson^^ that it is easy to de- 
termine whether a bar (bond) is redundant with respect 
to the bonds that are already in the lattice. 




FIG. 3. The stressed (AS) backbone on a lattice of size 
L = 1024 with open boundaries. 

If the bonds of the lattice are replaced by Hooke's 
springs, a redundant bond leads to internal stresses in 
the spring lattice. A redundant bond causes a violation 
of the condition 3nbod—b = 3 by adding an extra spring to 
the lattice. The algorithm checks to see if 3ni,od — 6 = 3 is 
satisfied by using "bipartite matching" or "the pebble 
game" ^^'^^ to see if each of the bodies' degrees of freedom 
can be "matched" to the bonds of the configuration (note 
that JT use the "pebble game" in the original bar-joint 
representation) . 

If an extra or redundant spring is added to a cluster that 
is already rigid (i.e. already satisfies 3ni,od — b = S), the 
matching algorithm identifies the bonds which become 
internally stressed when the extra spring is added. We 
then give these internally stressed bonds the same cluster 
label. In this way we find "internally stressed (IS) clus- 
ters" . Finally, an applied tensile stress can be mimiced by 
adding two rigid beams to the two opposite sides of the 
lattice, and then by adding a fictitious bond between the 
rigid beams. In this way, we determine when the lattice 
can support an applied tensile stress. Full details of the 
algorithm will appear elsewhere^^. In JT, the bond prob- 
ability is fixed. Our algorithm is complimentary to theirs 
as we add bonds one at a time, so we find the percola- 
tion concentration exactly for each sample. We chose this 
method not only because it is very efficient (comparable 
to JT), but also so that we can find the exact backbone 
for each sample, and hence the "critical (or red) bonds" 
of the backbone (see below). 



There are several ways to define the onset of stress trans- 
mission through a lattice. The two which are most phys- 
ically appealing are: 

1. The point at which an applied stress (AS) is trans- 
mitted across the lattice and; 

2. The point at which internally stressed regions (IS) 
connect together to form stressed clusters of macro- 
scopic size. 

Both of these definitions have simple representations in 
terms of a lattice of Hooke's springs. The first (AS) cor- 
responds to a random Hooke's spring lattice to which, for 
example, a tensile stress is applied, while the second (IS) 
corresponds to the internal stresses in a random Hooke's 
spring lattice with random natural lengths. We study the 
stressed backbone of these lattices as a function of site di- 
lution. We also tested the effect of boundary conditions 
on these two definitions of rigidity percolation, because a 
local change in rigidity (e.g. by adding a bond) can be 
transmitted over long distances so boundaries might be 
more important in this problem than in connectivity per- 
colation. However, we find that in the large-lattice limit 
both the AS and IS percolation definitions lead to the 
same, boundary condition independent, threshold. This 
behavior is presented in Fig. 2, from which wc find that 
Pc — 0.6975±0.0003. On regular lattices, previous work^^ 
using direct solution of the force equations lead to esti- 
mates close to Pc = 0.713 for samples of up to size L = 75. 
As can be seen from Fig. 2, this is consistent with the 
result on random lattices, although at that lattice size, 
there is still considerable dependence on boundary con- 
ditions. However, in general there is no reason to believe 
that the percolation threshold on random lattices should 
be the same as that on regular lattices. This difference 
is illustrated by the configuration of Fig. lb. On a reg- 
ular lattice that configuration is not rigid to shear, but 
if the lattice sites are displaced, it becomes rigid. That 
is because on a regular lattice, the three bars are paral- 
lel, so these constraints are "degenerate". Thus for that 
configuration, the random lattice is more rigid than the 
regular lattice. 
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FIG. 4. a) Finite-size-scaling plot of the normalised num- 
ber of bonds Nb /L'^ on the stressed backbone at the percola- 
tion point for: AS with periodic b.c.'s (*); AS with open b.c.'s 
IS with periodic b.c.'s (X) and ; IS with open b.c.'s (O). 
b) Finite-size-scaling plot of the normalised number of "red" 
bonds on the backbone, Nr/LF' (AS with periodic b.c's). The 
lattice sizes used and number of configurations were the same 
as for Fig. 2. 

However it is easy to construct configurations which 
are more rigid on a regular lattice (e.g. a sequence of 
aligned bonds forming a "guy" wire), so it is unclear as 
to whether displaced lattices have a lower or higher pc 
than regular lattices. 

Prom the variation in the percolation concentration 
Apc ~ L~^^^ , we are able to find the correlation length 
exponent. We did this for three ways of defining Apc, 
namely: 

a) {pc{L) -pc(oo)); 

b) (p°P"" - pP'^"°dtc'^ and; 

c) v'(<Pc >- <Pc >^), 

and for several types of boundary conditions in each case. 
From these extensive calculations, we find v = 1.16±0.03. 
Although it is not the main focus of this paper, we note 
that in JT, it is claimed that the infinite cluster. Poo 
(which includes internally stressed bonds (stressed back- 
bone), and unstressed bonds which satisfy 2n — b = 3) 
has a fractal dimension around Df ~ 1.86. 
If we assume a second order behavior in Pqo , we find a 
similar fractal dimension. However, a mean field theory^ 
suggests that the rigidity transition is first order (so 
Df = 2 in 2-d), and similarly on Bethe lattices the rigid- 
ity transition is first order^^. Thus we have also tested 
the possibility of a first order transition^^ in P^, and 
find that the data is consistent with a weakly first-order 
transition, with the first-order jump APqo ~ 0.085 at Pc- 
However, even larger lattices (up to of order L = 10, 000) 
are needed to determine convincingly whether, in 2-d, Pqo 
is first order. 

An example of a stressed backbone at the percolation 
point is presented in Fig. 3. We measured the number 
of bonds on backbones such as that shown in Fig. 3, and 
the results of a scaling plot are presented in Fig. 4. From 
this figure, we find Db = 1.78 ± 0.02. This backbone 
dimension is different than that for connectivity percola- 
tion where the backbone dimension is 1.62 ± 0.01, and it 
is also considerably larger than that of the stressed back- 



bone of regular triangular lattices 1.64 ± 0.05 found by 
direct solution* on small lattices (up to L = 80). The 
latter discrepancy could be due to a fundamental differ- 
ence between the random and regular lattices, but it also 
could be due to imprecise estimates of the percolation 
concentration in previous work due to finite size effects 
(see Fig. 2). At the percolation point, there are a set 
of bonds whose removal leads to loss of backbone rigid- 
ity. We calculated the number of these critical red bonds, 
iV/j, and their scaling behavior at the percolation point 
is also shown in Fig. 4. They scale as Nr ~ L^, with 
X = 0.85±0.05. This is consistent with x = although 
we have no analytic argument for why this should be so. 
Previous work on the elastic exponents of regular trian- 
gular, central force, networks have produced conflicting 
results. Although the early work^ gave an exponent in 
the range, 1.3 < f/v < 2.0, later work suggested that 
the central force and bond-bending (angular force) prob- 
lems were in the same universality class*' so that"'^^ 
f /v ~ 3.0. There have even been suggestions^^'^^ that in 
2-d, site percolation has exponent near f/v^ 1.0, while 
bond percolation has exponent near j jv ^ 3.0. Finally, 
there is a recent mean field theory^ which gives exponent 
/~1.5. 

The difficulty in obtaining good estimates have been as- 
cribed to: a) unusually strong accumulation of roundoff 
errors'', and b) lack of precision* in the estimate of Pc- 
We find that roundoff errors are largely eliminated if we 
use the graph theory method to remove all non-stressed 
bonds before applying the conjugate gradient method. 
In addition we know pc exactly for each configuration, so 
we do not have to study a range of p using an interative 
solver. Thus we have been able to study the elastic con- 
stants for lattice sizes which were previously inaccessible 
(up to linear size L = 512 - see Fig. 5). 
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FIG. 5. Finite-size-scaling behavior of the Young's modu- 
lus for: regular lattices (+); lattices randomly displaced by 
0.2 (A) and; lattices with random bond angles (□). The 
lattice sizes used (number of configurations) were as follows; 
16(20,000), 32(10,000), 64(1,000 - 2,000), 128(100 - 300), 
256(20-40), 512(4-6). 

As expected, a certain number of the "generic" back- 
bones are not rigid on a regular lattice due to degenera- 
cies. However, the fraction of the backbones that are non- 
rigid on regular lattices increases very slowly with lattice 
size, and is still only 50% at L = 512. Now note that 



if the sites of a generic backbone, which is non-rigid on 
a regular lattice, are displaced by a small amount A, the 
elastic modulus of that backbone is O(A^). Thus, the 
elastic constants of generic backbones are usually non- 
universal, for sizes accessable to simulation, even for lat- 
tices displaced by 0.2 (see Fig. 5). 

To avoid the slow size effect caused by proximity to the 
regular lattice limit, we also studied a model where the 
locations of the elements of the elastic matrix were set by 
the connectivity of the backbone. To mimic the highly 
displaced lattice (large A) limit, we assign each bond an 
angle to the x-axis which is drawn from a random distri- 
bution of angles (on the interval [0, 360], and calculate the 
elastic constant using these angles in the force equations. 
The results for this "random angle model" are also shown 
in Fig. 5 (each present bond has unit spring constant). 
We found that the value f jv ^ 1.45 is quite universal in 
this limit. 
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